%%% programa para procesar los datos de consumo diario
clear all
tic
load codensa.mat%% contiene los datos en clientes1_650188 como una matriz de 600000x39, donde del
%%%campo 4 al 38 se encuentran los datos de interes. la columna 39 es el
%%%tipo donde 1 es residencial, 2 es comercial, 3 es industrial y 4 es
%%%oficial
datacons=codensa(:,4:38);%%datos de consumo
%%%Borrar Nans de datacons
for i=1:size(datacons,2)
    clear temp
    temp=isfinite(datacons(:,i));
datacons=datacons(temp,:);
end
%%% calculo de media y std
[A]=[mean(datacons,2) std(datacons,0,2)];
%%%% Borrar los Nan
A2=A(isfinite(A(:,1)),:);
%%% Trabajar con Varianza normalizada respecto al promedio
A3=(A2(:,2)./A2(:,1));
    A3=[A2 A3]; 
%%% Si el promedio es cero eliminar los datos Nan de A3
 A4=A3(find(A3(:,1)~=0),:);%%% columna 1 es la media, columna 2 es la varianza, col3 es la varianza relativa
 
 %%Juntar desviaciones y promedios con estructira nro cuenta,
 %%estrato,clase,mean,desvest,desvestrelativa
 
 A7=[codensa(find(A3(:,1)~=0),[1 2 39]) A4];
 %%%% Calcular el coeficiente de correlaci?n de media con desviacion
 %%%% estandar relativa por tipo de cliente
 
 %%%Corrleacion para Residenciales
 A_res=corrcoef(A7(find(A7(:,3)==1),4:6))
 
 %%Correlacion para Comerciales
A_com=corrcoef(A7(find(A7(:,3)==2),4:6))
 %%COrrelacion para industriales
 A_ind=corrcoef(A7(find(A7(:,3)==3),4:6))
 %%COrrelacion para oficiales
A_off=corrcoef(A7(find(A7(:,3)==4),4:6))

%%% Media y desviacion estandar por estrato para residenciales
for i=1:6
 
     medvarianzaestrato(i,:)=[i mean(A7(find(A7(:,3)==1 & A7(:,2)==i),4)) std(A7(find(A7(:,3)==1 & A7(:,2)==i),6))]; %%encuentra los que son residenciales y tienen consumo medio menor a 30diario o 30 mensual
 
end
%%% Comportamiento Varianza  por rangos CODENSA
%%%Rangos Residenciales 300,550,900
for i=1:4
 if i==1
     varianza_res(i,:)=[i mean(A7(find(A7(:,3)==1 & A7(:,4)<=300/30),6)) std(A7(find(A7(:,3)==1 & A7(:,4)<=300/30),6))]; %%encuentra los que son residenciales y tienen consumo medio menor a 30diario o 30 mensual
 elseif i==2
     varianza_res(i,:)=[i mean(A7(find(A7(:,3)==1 & A7(:,4)>300/30 & A7(:,4)<=550/30),6)) std(A7(find(A7(:,3)==1 & A7(:,4)>300/30 & A7(:,4)<=550/30),6))]; 
 elseif i==3
     varianza_res(i,:)=[i mean(A7(find(A7(:,3)==1 & A7(:,4)>550/30 & A7(:,4)<=900/30),6)) std(A7(find(A7(:,3)==1 & A7(:,4)>550/30 & A7(:,4)<=900/30),6))]; 
 elseif i==4
     varianza_res(i,:)=[i mean(A7(find(A7(:,3)==1 & A7(:,4)>900/30),6)) std(A7(find(A7(:,3)==1 & A7(:,4)>900/30),6))]; 
 end

end
bar(varianza_res(:,1),varianza_res(:,2));title('Comportamiento Varianza Usuarios Residenciales seg?n rangos CODENSA');xlabel('Rangos');ylabel('Kwh/d?a');


%%%Rangos Industriales 300,660,1100,5000,12000,50000,1000000
for i=1:7
 if i==1
     varianza_ind(i,:)=[i mean(A7(find(A7(:,3)==3 & A7(:,4)<=300/30),6)) std(A7(find(A7(:,3)==3 & A7(:,4)<=300/30),6))]; %%encuentra los que son residenciales y tienen consumo medio menor a 30diario o 30 mensual
 elseif i==2
     varianza_ind(i,:)=[i mean(A7(find(A7(:,3)==3 & A7(:,4)>300/30 & A7(:,4)<=660/30),6)) std(A7(find(A7(:,3)==3 & A7(:,4)>300/30 & A7(:,4)<=660/30),6))]; 
 elseif i==3
    varianza_ind(i,:)=[i mean(A7(find(A7(:,3)==3 & A7(:,4)>660/30 & A7(:,4)<=1100/30),6)) std(A7(find(A7(:,3)==3 & A7(:,4)>660/30 & A7(:,4)<=1100/30),6))]; 
 elseif i==4
     varianza_ind(i,:)=[i mean(A7(find(A7(:,3)==3 & A7(:,4)>1100/30 & A7(:,4)<=5000/30),6)) std(A7(find(A7(:,3)==3 & A7(:,4)>1100/30 & A7(:,4)<=5000/30),6))];
 elseif i==5
     varianza_ind(i,:)=[i mean(A7(find(A7(:,3)==3 & A7(:,4)>5000/30 & A7(:,4)<=12000/30),6)) std(A7(find(A7(:,3)==3 & A7(:,4)>5000/30 & A7(:,4)<=12000/30),6))];
 elseif i==6
     varianza_ind(i,:)=[i mean(A7(find(A7(:,3)==3 & A7(:,4)>12000/30 & A7(:,4)<=50000/30),6)) std(A7(find(A7(:,3)==3 & A7(:,4)>12000/30 & A7(:,4)<=50000/30),6))];
 elseif i==7
     varianza_ind(i,:)=[i mean(A7(find(A7(:,3)==3 & A7(:,4)>50000/30 & A7(:,4)<=1e6/30),6)) std(A7(find(A7(:,3)==3 & A7(:,4)>50e3/30 & A7(:,4)<=1e6/30),6))];
 end
end
figure
bar(varianza_ind(:,1),varianza_ind(:,2));title('Comportamiento Varianza Usuarios Industriales seg?n rangos CODENSA');xlabel('Rangos');ylabel('Kwh/d?a');

 
 %%% Calculo del VAR
 %%% Calculo de diferencias porcentuales
 deltarel=100*diff(datacons,1,2)./datacons(:,1:size(datacons,2)-1);
 %%%OJO SOLO PARA PRUEBAS!!!!!
 %deltarel=deltarel(1:10,:);
 %%%% Hay que tener cuidado con los datos de consumo Inf en deltarel ya que
 %%%% se producen cuando en un periodo hay consumo 0. Para evitar este
 %%%% problema se convierten todos los inf de deltarel en NaN para que as?
 %%%% la funci?n percentil no los tome en cuenta
 deltarel(find(deltarel==Inf))=NaN;
 
 VAR=[prctile(deltarel,5,2) prctile(deltarel,95,2)];
 umbrales=[codensa(temp,1:3) VAR];
toc
 warning('Listo el calculo del VAR')
 tic
 
 %%%Prueba utilizando criterios de VAR para detecci?n de anomal?as
 %% Esta prueba se realizar? utilizando los valores de los ultimos ke meses
 %% para verificar la efectividad del metodo
 %%deltarel tiene longitud de 34 campos
 tic
 ke=12;%%% periodo de evaluaci?n ultimos ke meses
 %%Matriz de sospechosos presumiendo inocencia con valores 0
 sospechosossinventana=zeros(length(deltarel),ke);%%% 
 for i=1:ke
     %%%Calcular los umbrales
     kep=size(deltarel,2)-ke-1+i;
     VARi=[prctile(deltarel(:,1:kep),5,2) prctile(deltarel(:,1:kep),95,2)];
     %%%Verificar si los nuevos incrementos superan los umbrales comparando
     %%%deltarel con Vari
%      for j=1:length(deltarel)
%          if deltarel(j,kep+1)<VARi(j,1) | deltarel(j,kep+1)>VARi(j,2)%% si supera los umbrales
%              sospechosos(j,i)=1;
%          else
%              sospechosos(j,i)=0;
%          end
%          
%          
%      end

%%%Comparara los nuevos valores con los limites inferiores 5% y
%%%superioireres 95% mediante una diferencia

complimite5=deltarel(:,kep+1)-VARi(:,1);%%comparacion con VAR5%
complimite95=deltarel(:,kep+1)-VARi(:,2);%%comparacion con VAR5%
%%% Sospechosos
sospechosossinventana(find(complimite5<=0),i)=1;%%tomar como sospechos cuando es menor que VAR5
sospechosossinventana(find(complimite95>=0),i)=1;%%tomar como sospechosos cuando es mayor a VAR95

     i
 end
 %%%Organizar clientes y sospechosos
 criticavar=[codensa(temp,[1 2 39]) sospechosossinventana];  
 toc
 
%%%Calculo de Var con ventanas
%%%Ventana de 12 meses
vent=12;
 sospechososvent12=zeros(length(deltarel),ke);%%% 
 for i=1:ke
     %%%Calcular los umbrales
     kep=size(deltarel,2)-ke-1+i;
     VARi=[prctile(deltarel(:,kep-vent+1:kep),5,2) prctile(deltarel(:,kep-vent+1:kep),95,2)];
     %%%Verificar si los nuevos incrementos superan los umbrales comparando
     %%%deltarel con Vari
%      for j=1:length(deltarel)
%          if deltarel(j,kep+1)<VARi(j,1) | deltarel(j,kep+1)>VARi(j,2)%% si supera los umbrales
%              sospechosos(j,i)=1;
%          else
%              sospechosos(j,i)=0;
%          end
%          
%          
%      end

%%%Comparara los nuevos valores con los limites inferiores 5% y
%%%superioireres 95% mediante una diferencia

complimite5=deltarel(:,kep+1)-VARi(:,1);%%comparacion con VAR5%
complimite95=deltarel(:,kep+1)-VARi(:,2);%%comparacion con VAR5%
%%% Sospechosos
sospechososvent12(find(complimite5<=0),i)=1;%%tomar como sospechos cuando es menor que VAR5
sospechososvent12(find(complimite95>=0),i)=1;%%tomar como sospechosos cuando es mayor a VAR95

     i
 end
 %%%Organizar clientes y sospechosos
 criticavar12=[codensa(temp,[1 2 39]) sospechososvent12];  
 toc
 
 
 
%%%Calculo de Var con ventanas
%%%Ventana de 18 meses
vent=18;
 sospechososvent18=zeros(length(deltarel),ke);%%% 
 for i=1:ke
     %%%Calcular los umbrales
     kep=size(deltarel,2)-ke-1+i;
     %deltarel(:,kep-vent+1:kep)
     VARi=[prctile(deltarel(:,kep-vent+1:kep),5,2) prctile(deltarel(:,kep-vent+1:kep),95,2)];
     %%%Verificar si los nuevos incrementos superan los umbrales comparando
     %%%deltarel con Vari
%      for j=1:length(deltarel)
%          if deltarel(j,kep+1)<VARi(j,1) | deltarel(j,kep+1)>VARi(j,2)%% si supera los umbrales
%              sospechosos(j,i)=1;
%          else
%              sospechosos(j,i)=0;
%          end
%          
%          
%      end

%%%Comparara los nuevos valores con los limites inferiores 5% y
%%%superioireres 95% mediante una diferencia

complimite5=deltarel(:,kep+1)-VARi(:,1);%%comparacion con VAR5%
complimite95=deltarel(:,kep+1)-VARi(:,2);%%comparacion con VAR5%
%%% Sospechosos
sospechososvent18(find(complimite5<=0),i)=1;%%tomar como sospechos cuando es menor que VAR5
sospechososvent18(find(complimite95>=0),i)=1;%%tomar como sospechosos cuando es mayor a VAR95

     i
 end
 %%%Organizar clientes y sospechosos
 criticavar18=[codensa(temp,[1 2 39]) sospechososvent18];  
 toc
 
 
%%%Calculo de Var con ventanas
%%%Ventana de 6 meses
vent=6;
 sospechososvent6=zeros(length(deltarel),ke);%%% 
 for i=1:ke
     %%%Calcular los umbrales
     kep=size(deltarel,2)-ke-1+i;
     VARi=[prctile(deltarel(:,kep-vent+1:kep),5,2) prctile(deltarel(:,kep-vent+1:kep),95,2)];
     %%%Verificar si los nuevos incrementos superan los umbrales comparando
     %%%deltarel con Vari
%      for j=1:length(deltarel)
%          if deltarel(j,kep+1)<VARi(j,1) | deltarel(j,kep+1)>VARi(j,2)%% si supera los umbrales
%              sospechosos(j,i)=1;
%          else
%              sospechosos(j,i)=0;
%          end
%          
%          
%      end

%%%Comparara los nuevos valores con los limites inferiores 5% y
%%%superioireres 95% mediante una diferencia

complimite5=deltarel(:,kep+1)-VARi(:,1);%%comparacion con VAR5%
complimite95=deltarel(:,kep+1)-VARi(:,2);%%comparacion con VAR5%
%%% Sospechosos
sospechososvent6(find(complimite5<=0),i)=1;%%tomar como sospechos cuando es menor que VAR5
sospechososvent6(find(complimite95>=0),i)=1;%%tomar como sospechosos cuando es mayor a VAR95

     i
 end
 %%%Organizar clientes y sospechosos
 criticavar6=[codensa(temp,[1 2 39]) sospechososvent6];  
 toc
 
 
% %%%Calculo de Var con ventanas
% %%%Ventana de 24 meses
% vent=24;
%  sospechosos24=zeros(length(deltarel),ke);%%% 
%  for i=1:ke
%      %%%Calcular los umbrales
%      kep=size(deltarel,2)-ke-1+i;
%      VARi=[prctile(deltarel(:,kep-vent:kep),5,2) prctile(deltarel(:,kep-vent:kep),95,2)];
%      %%%Verificar si los nuevos incrementos superan los umbrales comparando
%      %%%deltarel con Vari
% %      for j=1:length(deltarel)
% %          if deltarel(j,kep+1)<VARi(j,1) | deltarel(j,kep+1)>VARi(j,2)%% si supera los umbrales
% %              sospechosos(j,i)=1;
% %          else
% %              sospechosos(j,i)=0;
% %          end
% %          
% %          
% %      end
% 
% %%%Comparara los nuevos valores con los limites inferiores 5% y
% %%%superioireres 95% mediante una diferencia
% 
% complimite5=deltarel(:,kep+1)-VARi(:,1);%%comparacion con VAR5%
% complimite95=deltarel(:,kep+1)-VARi(:,2);%%comparacion con VAR5%
% %%% Sospechosos
% sospechososvent24(find(complimite5<=0),i)=1;%%tomar como sospechos cuando es menor que VAR5
% sospechososvent24(find(complimite95>=0),i)=1;%%tomar como sospechosos cuando es mayor a VAR95
% 
%      i
%  end
 %%%Organizar clientes y sospechosos
% criticavar24=[codensa(temp,[1 2 39]) sospechososvent24];  
 % toc
 
 
%  %%% Calculo de Kolmogorv smirnov a valores extremos
%  for i = 1:size(datacons,1)
%      evpar=evfit(datacons(i,:));%%Calcula los parametros de valores extremos
%      if evpar(1)==0 & evpar(2)==0 %%% Si los parametros son cero kolmogorov es 0
%          kolm(i)=0;
%      else  
%      evdist=evcdf(datacons(i,:),evpar(1),evpar(2));%%calcula la distribucion acumulada de valores extremos
%      kolm(i)=kstest(datacons(i,:),[datacons(i,:)' evdist']);
%      kolm_norm(i)=kstest(datacons(i,:));
%      end
%      clear evpar
%      clear evdist
%  end
%  toc
 %save varcodensasim

save varcodensaventanas
